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We present first-principles calculations of the structural and electronic properties of Si(001)-SiO2 
interfaces. We first arrive at reasonable structures for the c-Si/a-Si02 interface via a Monte-Carlo 
simulated annealing applied to an empirical interatomic potential, and then relax these structures 
using first-principles calculations within the framework of density-functional theory. We find a tran- 
sition region at the interface, having a thickness on the order of 20 A, in which there is some oxygen 
deficiency and a corresponding presence of sub-oxide Si species (mostly Si +2 and Si +3 ). Distribu- 
tions of bond lengths and bond angles, and the nature of the electronic states at the interface, are 
investigated and discussed. The behavior of atomic oxygen in a-SiC>2 is also investigated. The per- 
oxyl linkage configuration is found to be lower in energy than interstitial or threefold configurations. 
Based on these results, we suggest a possible mechanism for oxygen diffusion in a-Si02 that may be 
relevant to the oxidation process. 
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I. INTRODUCTION 

Understanding the atomic structure of ultrathin sili- 
con dioxide films on Si(100) substrates is an outstand- 
ing problem-||Of great importance for microelectronic 
applications.lira The thickness of these films in commer- 
cial devices has dropped well below 100A. Film thick- 
nesses of 40A and below are now being explored in ex- 
perimental devices. Although such ultrathin oxide films 
are now in wide use, and considerable information about 
them has been amassed-using a variety of different exper- 
imental techniques ]3un there is little consensus about 
the microscopic structure and oxidation kinetics of the 
films. Among the factors that have so far hindered de- 
tailed studies of the interface are the amorphous nature 
of the Si02 region, the difficulty of probing a buried inter- 
face, and the dependence of sample properties on prepa- 
ration conditions. As a result, several issues relating to 
the detailed chemical structure of the interface, and to 
the oxidation mechanism of thin oxide films, continue to 
be controversial. 

Photoemission spectroscopy (PES) and photoelectron 
spectroscopy (XPS) have been the major experimental 
techniques used to study the interface structure. These 
experiments indicate the presence of a transition region 
of sub-oxidized Si near the interface, with all three in- 
termediate partial oxidation states observed. However, 
different experiments suggest conflicting pictures of the 
structure. For example, PES experiments .indicate that 
the total suboxide is about two monolayers ,lH while XPS 
results suggest it is one monolayer thick. □ Using medium- 
energy ion-scattering spectroscopy (MEIS) to study the 
interface, Gusev et a/Ej and Lu et al&3 have studied se- 
quential isotopic exposures of oxygen and found that the 
growth mechanism of interfaces depends on the thick- 
ness of the films. Further complicating the task has been 



the lack of suitable interface models which serve as ref- 
erence points for analysis of spectr oscop ic results. The- 
oretical works by Pasquarello et alh3t3 have attempted 
to construct interface models by attaching tridymite, a 
crystalline form of SiC>2, to Si(100). Although different 
suboxide states can be created, one would expect such an 
artificial model would introduce significant stress across 
the interface. A recent— experimental work suggest an 
abrupt interface model,Ej as opposed to the graded in- 
terface suggested from earlier work.Ercl Additional uncer- 
tainty has been generated by some other studies *t|f*|t 
the interpretation of the sub-oxide peaks observed. 

Another aspect that is of great interest is the oxida- 
tion kinetics for thin films. The formation of thick oxide 
films (> 200 A) is well -described phenomenologically by 
the Deal-Grove modeljill in which molecular oxygen dif- 
fuses to the SiO^/Si interface and reacts with silicon at 
the interface. This model predicts a linear relationship 
between oxide thickness and oxidation time. However, it 
is well known that the Deal-Grove model breaks down for 
the case of ultrathin films (<100A); in this case, the oxi- 
dation kinetics have been shown to be fastep-jthaB, would 
be expected from the linear relationship.™ tttej How- 
ever, there is little consensus about the mechanisms of 
oxidation kinetics in this case— .Seyeral phenomenological 
models have been proposed .Ej til Some of them fit the 
experimental data on oxidation kinetics quite satisfacto- 
rily with a large number of fitting parameters. However, 
most of the models arc without direct experimental sup- 
port. An analysis of kinetic results alone does not allow 
one to distinguish among different models. A few key 
problems associated with the oxidation mechanism still 
remain unresolved. For example, it is uncertain exactly 
where oxidation takes place, and the mechanism of oxy- 
gen diffusion through a-Si02 region in not clear. 

In this paper, we choose to address two main issues, 
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namely, the microscopic structure of Si/Si02 near the 
interface and the behavior of atomic oxygen in a-SiC>2. 
Firstly, we arrive at plausible structures for the Si/SiC>2 
interface by performing Monte Carlo simulated anneal- 
ing using an empirical potential. Then we relax these 
structures and study their properties via first-principles 
density-functional calculations. We obtain the statistical 
distribution of sub-oxide species near the interface region, 
and find a graded transition region with a width of about 
10 A. Secondly, for excess atomic oxygen in a-SiC>2, we 
compare the energy of the interstitial and peroxyl bridge 
configurations. It is found that the latter is always lower 
in energy. This suggests a possible mechanism of oxygen 
hopping among neighboring Si-Si bonds. 

The plan of the paper is as follows. Section II gives 
a brief description of the technique used to perform the 
calculations. In Sec. Ill we present our work on the mi- 
croscopic structure of the Si/Si02 interface, and Sec. IV 
contains the results on oxidation kinetics of atomic oxy- 
gen in a-Si02- We summarize the work in Sec. V. 



II. METHODS 

Our theoretical analysis consists mainly of two parts. 
In the first part, we employ Metropolis Monte Carlo 
(MC) simulations together with an empirical model po- 
tential representing the structural energies to arrive, via 
simulated annealing, at candidate structures for the c- 
Si/a-Si02 interface. In the second part of the analysis, 
these candidate structures are relaxed and analyzed using 
first-principles plane-wave pseudopotential calculations. 
Because the structures of interest-involve coordination 
changes, we have followed HamannE3 in employing a gen- 
eralized gradient approximation (GGA) to the exchange- 
correlation potential, specifically that of Perdew, Burke 
and Ernzcrhof (PBE)EI 

For use in generating the c-Si/a-Si02 interface struc- 
tures within the MC approach, we designed a simple 
empirical-potential model for the structural energy that 
is roughly based on the formalism of Keating. Ell The only 
degrees of freedom that appear in this model are the Si 
atom coordinates; a direct Si-Si bond is designated as a 
"short" (S) bond, while a Si-O-Si bridge configuration is 
designated as a "long" (L) bond. In this way, the whole 
c-Si/a-Si02 structure is replaced by a network of "short" 
and "long" bonds between Si atoms. Its energy is taken 
to be 
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(1) 



The first term represents the bond stretching term, where 
and d$ take values kf = 9.08eV/A 2 and = 2.35A 
for a "short" bond, or k^ = 1.89eV/A 2 and = 3.04A 



for a "long" bond. Similarly, the second term repre- 

(ii) 

sents the bond-bending term, where kg takes values 
kl s = 3.58eV, /cf L = 3.81eV, and fc^ L = 4.03eV, where 
i,j are summed over neighboring bonds only. Here Cjj is 
the cosine of the angle subtended by bonds i and j, and 
cj\^ = cos(109.47°) = —1/3 is the corresponding refer- 
ence value. The parameters for the empirical potential 
were chosen such that they reproduce not only the equi- 
librium structures, but also the bulk and shear moduli, 
of crystalline Si and /3-cristobalite structures. The only 
parameter that cannot be fixed empirically is the bond- 
bending term kg L associated with a mixture of S and L 
bonds. However, it is found that kg S and k^ h are very 
close to each other, so kg h was chosen as the average of 
these. 

There are two types of discrete moves allowed in 
our MC simulations, namely, bond switching and bond 
conversion moves. For the bond switching moves, we 
adopted tlie algorithm of Wooten, Winer, and Weaire 
(WWW)a To randomize the structure, one of the sec- 
ond neighbors of a randomly chosen atom is switched 
to become its first neighbor, and that first neighbor be- 
comes the second neighbor (see Fig. [l]). Notice that the 
bonds being switched can be either S or L bonds. For 
the bond conversion moves, a randomly chosen S or L 
bond is swapped with a neighboring bond of opposite 
type. In all of the discrete moves, the total number of 
S and L bonds are conserved, meaning that the number 
of oxygen atoms in the system does not change during 
the MC simulations. After each move, the structure is 
relaxed within the Keating-like model, and the energy 
difference between the relaxed structures before and af- 
ter the move is computed. The Metropolis algorithm is 
then employed to decide whether the new structure is to 
be accepted. In the case of the Si/Si02 interface, certain 
atoms in the c-Si region are fixed (i.e., not allowed to re- 
lax, and no bond switching or conversion moves involving 
them are allowed). 

After generating structures from the MC simulations, 
first-principles plane-wave pseudopotential calculations 
are carried out to produce accurate and relaxed struc- 
tures. The exchange and correlation energy aad potential 
are implemented in the PBE-GGA scheme.E^I The pseu- 
dopotential for O is that used in Ref. 26, and was, gen- 



erated using an ultrasoft pseudopotential schemeO The 
ultrasoft pseudopotential for Si is generated with a core 
radius of 1.3 a.u. and valence states of 3s and 3p. Self- 
consistent total-energy and force calculations are used to 
relax the atomic coordinates until the forces are less than 
0.25 eV/A. A cutoff energy of 25 Ry is used for the calcu- 
lations in Sec. IV, but this was reduced to 15 Ry for the 
majority of total-energy calculations of interface struc- 
tures presented in Sec. III. The low energy cutoff for the 
interface structures was chosen so as to reduce the com- 
putational burden. It is found that when we raise the 
cutoff energy to 20 Ry, there is only a slight effect on 
the structural configurations. More specifically, there is 
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FIG. 1. The WWW bond-switching move (Ref.||) used in 
the Monte Carlo simulations, (a) Initial configuration. Atoms 
1 and 4 are first neighbors, and atoms 1 and 7 are second 
neighbors, (b) After bond switch. Atoms 1 and 7 are now 
first neighbors, while atoms 1 and 4 are second neighbors. 

less than 1% of error in bond lengths and 5% of error in 
bond angles when we compare structures generated with 
different cutoff energies. So, a cutoff energy of 15 Ry is 
enough to provide accurate results for our analysis. For 
all cases, the Brillouin zone sampling consists of a single 
k-point at T. This sampling choice is already represent- 
ing a high k-point density in the reciprocal space, as large 
periodic supercells are being used in the calculations. 

In order to treat the PBE-GGA exchange and corre- 
lation potential efficiently during the plane-wave calcu=. 
lations, we employ the approach of White and Bird.Ej 
In this method, by writing the energy and potential in a 
discretized form, there is no need to represent the charge 
density and its gradient on a finer real-space mesh. This 
not only reduces the computational time and memory re- 
quirements, but also improves the convergence and sta- 
bility of the calculations. 

III. THE Si/SiOa INTERFACE 

The initial input structure for the MC annealing pro- 
cedure is a crystalline supercell containing either 12 or 
40 layers of Si atoms stacked along the [001] direction, 
each layer consisting of 4 atoms in a cell of dimensions 
y/2a x \[2a in the x-y plane. The three central layers of 
Si atoms are held fixed in each structure, to form the por- 
tion that will remain as the c-Si region. The bonds in that 
region, plus either one or two layers of adjacent bonds, 
are initially assigned as S bonds; the remainder are as- 
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FIG. 2. Oxidation number of Si vs. distance of Si from 
interface, as obtained from the Monte Carlo simulations. 
Curves represent distribution profiles for five different oxi- 
dation states, smoothed by a gaussian broadening of width 
~2A. Data is compiled from six different structures with a 
total of ~10 3 Si atoms. 

signed as L bonds to form the "oxide." (The structure 
formed in this way is initially highly strained, but the ac- 
tion of the MC bond-switching moves, combined with the 
fact that the c-axis is allowed to relax on each MC step, 
soon anneals away most of this strain.) All bonds except 
those in the central three-layer c-Si region are allowed 
to be switched and/or converted during MC simulations. 
The MC annealing schedule starts with a temperature on 
the order of roughly 10 4 K with five sweeps. (Each sweep 
is a number of trial steps equal to the total number of 
atoms.) Then, the temperature is reduced by approx- 
imately a factor of two while the number of sweeps is 
increased by a factor of two for each consecutive simula- 
tion stage. The annealing schedule stops when the tem- 
perature reaches around 700K and the number of sweeps 
reaches the order of 10 2 . Further reducing the tempera- 
ture does not produce a significant effect in the resulting 
structures. 

In order to investigate the transition region near the 
interface as generated by the MC procedure, we use ini- 
tial input geometries having 40 layers of Si atoms, having 
a total of 160 Si atoms and 256 O atoms. The initial di- 
mensions of these periodic supercells are 7. 68 A x 7.68A x 
54.30A, but the third cell dimension relaxes to ~90A dur- 
ing the simulations. The resulting statistical distribution 
of Si oxidation number (i.e., the number of oxygen neigh- 
bors to a given Si atom) vs. distance from the interface 
is shown in Fig. |^. It shows that the size of the transi- 
tion region is about 20 JL._which agrees with some for- 
mer experimental work. EH The figure also shows that 
some Si +1 species are concentrated within ~5A of the 
interface, while some Si +2 species extend ~10A into the 
a-Si02 region. These cesults agree very well with some 
XPS and PES studies^ that find that Si +1 and Si+ 2 
states are localized within 6-10A of the interface. These 
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FIG. 3. Two examples of relaxed interface structures, 
as generated by Monte Carlo annealing followed by 
first-principles relaxation. Periodic replicas arise because the 
field of view spans two supercell dimensions in each of the 
directions parallel to the interface. 

experiments also suggest that Si +3 species extend ~30A 
into the bulk SiC>2, again in qualitative agreement with 
our simulations. Himpsel et. ala have proposed that 
these characteristic protrusions of Si +3 into the bulk SiC>2 
could be seen as the cores of misfit dislocations, which are 
caused by the large lattice mismatch of Si and SiC>2. 

Of course, we emphasize that our statistical results re- 
ported in Fig. |^ are generated purely from the MC sim- 
ulations, with no input from the first-principles calcula- 
tions. The z-positions of the atoms would shift slightly 
when relaxed with the first-principles calculations, but 
this typically shifts the points in the figure by less than 
0.2A. More significantly, one should ideally recompute 
the statistics using the first-principles total energies in 
the MC procedure, but such an approach would be pro- 
hibitively expensive. Nevertheless, the qualitative agree- 
ment with experiment regarding the distribution of sub- 
oxide species suggests that the interface structures as 
generated by the empirical MC procedure are fairly real- 
istic. 

First-principles calculations were only carried out for 
the structures with 12 layers of Si atoms (a total of 96 or 
112 atoms, depending on how many layers were oxidized). 
For the 40-layer structures, such calculations would be 
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FIG. 4. Histograms of Si-O-Si bond angles, for bond an- 
gles located (a) within 5 A of the interface, and (b) between 
5 and 10 A from the interface. 



too computationally demanding. The initial input di- 
mensions of the periodic supercell are 7.68Ax 7.68A x 
16.29A, but this relaxes to 7.68A x 7.68A x ~23A when 
the third dimension is allowed to relax during the MC 
procedure. First-principles calculations are then carried 
out on seven structures to relax both the atomic coordi- 
nates and the supercell z dimension, the latter typically 
being found to change by less than 6% from the empirical- 
potential MC value. Note that these seven structures 
were generated using the same MC procedure, but with 
different random number seeds, so that they can be re- 
garded as a small subset of the full ensemble that would 
be generated by the MC approach. 

Two examples of the relaxed interface structures are 
shown in Fig. [| Because of the lattice mismatch, a re- 
gion of pure a-Si02 cannot easily be attached directly to 
a c-Si substrate. As a result, as illustrated in the ex- 
amples of Fig. |^, sub-oxide species are numerous in the 
region near the interface. One might then expect that 
the structure in the oxide should differ significantly from 
that of bulk a-Si02 in this transition region. Indeed, 
our results support this. For example, Fig. |I] shows the 
statistical distribution of oxygen bond angles in differ- 
ent regions away from the interface. These distributions 
become broader and more evenly distributed as one gets 
further from the interface, eventually resembling fairly 
closely the accepted bulk a-Si02 distribution^ 

Another structural aspect we have examined is the Si- 
Si bond length in the Si02 region. We find that the 
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bond lengths become shorter when the average oxidation 
state of the two participating Si atoms increases from 
zero to two. The bond length is about 4% shorter than 
that of c-Si (2.35A) when the average oxidation is two. 
As the average oxidation state of the two Si increases 
to three, the bond length rebounds slightly, but is still 
about 3% shorter than that of the c-Si. Since an Si-0 
bond has less covalent character than a Si-Si bond, the 
initial shortening of the bond can be understood as a 
strengthening of the covalent character of the remaining 
Si-Si bonds, while the rebound for Si +3 is presumably 
related to the Coulomb repulsion between the increas- 
ingly positively charged Si atoms as they lose electrons to 
their more electronegative O neighbors. The Si-0 bond 
length decreases from about 1.67A to about 1.62A when 
the oxidation state of Si involved increases from one to 
four. As expected, the Si +4 -0 bond length is very close 
to the observed value (1.6A) in a-Si02- 

We have also investigated the electronic states of these 
interface structures. We do not find any defect state in 
the energy gap of our structures. This is not surprising, 
and reflects the fact that our geometries correspond to a 
continuous random network structure without any bro- 
ken bonds being present. We also specifically inspected 
the highest occupied molecular orbitals (HOMOs) and 
lowest unoccupied molecular orbitals (LUMOs) of these 
structures, but did not find any simple pattern for the 
localization of these orbitals. Sometimes the HOMO or 
LUMO states are confined mainly to the c-Si region (typ- 
ically being concentrated either at the top or bottom in- 
terface but remaining relatively delocalized in the x-y 
plane), but they also frequently have substantial ampli- 
tude at certain localized sites in the oxide region. We 
tried to determine whether abnormal bond lengths or 
bond angles might be associated with the localization of 
the LUMO or HOMO orbitals, but we found no simple 
correlation of this kind. 

In summary, in this part of our work we have shown 
that it is possible to generate reasonable interface struc- 
tures for the Si/Si02 interface despite the huge lattice 
mismatch, provided that one allows for the presence of 
sub-oxide states of Si in the transition region. The dis- 
tribution of sub-oxide states that emerges from our sim- 
ulations is generally consistent with experiments. More- 
over, our interface structures are entirely free of dangling 
bonds. Of course, by its construction the MC proce- 
dure automatically generates continuous random network 
(CRN) structures that are free of dangling bonds, but 
we then relaxed these structures using the first-principles 
calculations. While bonds might have broken or reformed 
during this relaxation procedure, such events were never 
observed in our calculations. That is, the first-principles 
relaxations always preserved the identity of the MC- 
generated CRN structures. This result is not trivial, and 
provides strong support for the robustness of the MC 
procedure. The first-principles calculations also provided 
Si-0 bond-length and Si-O-Si bond-angle distributions 
that are consistent with the the accepted picture for bulk 




FIG. 5. Possible bonding configurations for an extra oxy- 
gen atom in a-SiC>2. (a) "Threefold" configuration in which 
an extra oxygen binds to a bridge-bonded oxygen, (b) "Per- 
oxyl bridge" configuration in which an extra oxygen inserts 
into a bridge bond. Filled and open circles represent oxygen 
and silicon atoms respectively. 

a-Si02, and the electronic spectrum was checked for the 
absence of gap states. For all these reasons, we argue 
that the two-step procedure proposed here (MC genera- 
tion followed by first-principles relaxation) produces re- 
alistic interface structures. 



IV. OXIDATION KINETICS IN a-Si0 2 

In this part of the calculations, the energetics of atomic 
oxygen in an a-Si02 environment is investigated. Thus, 
for the MC simulations, only a random network of "long" 
bonds is needed as input. Periodic supercells (lO.lOAx 
lO.lOAx 7.14A) containing 49 atoms (including the ex- 
tra interstitial O atom) are used. The extra interstitial 
oxygen atom is first placed randomly inside the cavity 
formed by the a-si02 network, and the system is allowed 
to relax. We find that the extra oxygen atom binds to a 
neighboring oxygen atom that is already part of a Si-O- 
Si bridge bond, as illustrated in Fig. 5(a). The distance 
between the neighboring oxygens is only ^1.54 to 1.59 A, 
so that the bridge-bonded oxygen is essentially threefold 
coordinated. 

This threefold structure is a local minimum of the en- 
ergy landscape, but to determine whether it is truly sta- 
ble or only metastable, we then tried rearranging the two 
oxygen atoms into a peroxyl bridge structure as shown 
in Fig. 5(b). Again the structure is allowed to relax, 
and a local energy minimum associated with the peroxyl 
bridge configuration is found, with the distance between 
the two oxygen atoms being ~1.50 to 1.57A. The exis- 
tence of peroxyl bridge configurations in Si02 has already, 
been suggested on the basis of some experimental workE2l 
in which UV absorption evidence for peroxyl bridges was 
found in the presence of excess oxygen. Comparing the 
energies of the threefold and peroxyl configurations, we 
find that the peroxyl bridge is always the more stable 
configuration, with the energy difference varying in the 
range of about 0.9 to 2.8eV. We also find that there is 
a strong correlation between the relative stability of the 
peroxyl configuration and the length of the relaxed Si-Si 
distance in the peroxyl bridge; thus, the peroxyl is (not 
surprisingly) most likely to form where there is room for 
the two Si atoms to separate to make room for it. 
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FIG. 6. Proposed mechanism for diffusion of atomic oxygen 
in a-SiC>2. (a) Initial configuration with peroxyl bridge at left, 
(b) Peroxyl bridge excites into a threefold structure, (c) Extra 
oxygen atom jumps to a neighboring bond, (d) Threefold 
structure decays to form new peroxyl bridge at right. 

We find almost no energy barrier between the threefold 
and peroxyl configurations. This suggests that the three- 
fold configuration is only barely metastable (i.e., almost 
unstable). There could be two implications to be drawn 
from this result. First, the peroxyl bridge configuration 
can assist in the diffusion kinetics of oxygen. Sponta- 
neous formation and disintegration of peroxyl structures 
can occur during the oxidation process, and thus helping 
in promoting the diffusion of oxygen. Second, the per- 
oxyl bridge structure can hop among neighboring bonds 
of Si. Starting from the peroxyl bridge configuration, the 
hopping could occur by an initial formation of a threefold 
configuration, then a jump of the extra oxygen atom to a 
neighboring bond, and then the formation of a new per- 
oxyl bridge in that bond. This proposed hopping mech- 
anism is illustrated in Fig. |^. The precise energy barrier 
for this mechanism is not computed, as it depends in de- 
tail on the specific environment. However, we expect it 
should be of the same order of magnitude as the differ- 
ence of the peroxyl bridge and threefold configurations. 

This hopping mechanism is similar in spirit to one 
proposed recently by HamannJ22l who studied a per- 
oxyl bridge structure within an a-quartz supercell. In 
Hamann's model, the overall hopping process is similar 
to ours, except that in place of the intermediate structure 
that we propose in Figs. ^|(b) and ||(c), he has proposed 
an intermediate structure in which the extra oxygen atom 
is attached to the shared Si atom and acts as its fifth 
neighbor. The process then becomes a three-step mech- 
anism. It begins with a peroxyl bridge, then evolves into 
a fivefold Si configuration, and finally becomes a new 
peroxyl bridge in another bond. Our energy difference 
between the consecutive structures in the hopping mech- 
anism is roughly similar to that of Hamann's. 



V. SUMMARY 

We have studied a number of supercells to model 
the Si(100)/a-SiC>2 interface. A Monte Carlo simulated- 
annealing method, applied to a simple Keating-like model 
potential in which "long" and "short" bonds represent 
Si-O-Si and Si-Si bonds, is used to generate reasonable 
structural models of the interface. These structures are 
consistent with experimental observations of sub-oxide 
distributions. In particular, Si +1 and Si +2 sub-oxide 
states are concentrated within about 10A of the inter- 
face, with the Si +3 extends up to about 20A into the 
bulk SiC>2. As expected, the SiC>2 region shows more 
of the bulk a-Si02 character as one moves away from 
the interface. We have also investigated the HOMO and 
LUMO states of the interface structures. We confirmed 
that no deep gap states exist, but did not find any con- 
sistent pattern for the nature of the HOMO and LUMO 
orbitals. 

The kinetics of atomic oxygen in a-Si02 were also in- 
vestigated. The peroxyl bridge structure was found to be 
stable, with its energy computed to be lower than that 
of interstitial of threefold configurations. We suggest a 
possible mechanism for transport of atomic oxygen in 
the a-Si02 film, involving diffusion of the peroxyl bridge 
configuration, that could play a role in the interface oxi- 
dation during the formation of thin films. 
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